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We discuss a mechanical model which mimics the main features of the radiation matter interaction 
in the black body problem. The pure classical dynamical evolution, with a simple discretization of 
the action variables, leads to the Stefan- Boltzmann law and to the Planck distribution without any 
additional statistical assumption. 



The problem of black body radiation occupies a cen- 
tral role in the history of physics. The long debate 
around its properties stimulated a profound revision of 
old and firmly established concepts and eventually led to 
the birth of quantum mechanics. Yet it is somehow un- 
fortunate that such a profound revision of a fundamental 
theory has been based on a problem of so great complex- 
ity. Indeed, from the dynamical viewpoint, one is dealing 
with a nonlinear system with infinite degrees of freedom. 
Needless to say, such an infinite set of nonlinear differen- 
tial equations is not solvable. What is more important, 
and perhaps less noticed, is that even the statistical de- 
scription do not rest on solid grounds. As a matter of 
facts, classical ergodic theory which is now quite well de- 
veloped, is valid for systems with a finite number N of 
degrees of freedoms. A main difficulty stems from the 
fact that the two limits N — > oo and t — - > oo do not 
commute. The quantum statistical description is even 
more complex since it involves the additional difficulty 
that also the two limits t — > oo, h — > do not commute. 

In view of the relevance of the black body problem, 
it is perhaps worthwhile, after hundred years, a re- 
examination of the statistical properties with the help of 
modern computers and at the light of the recent progress 
in the study of nonlinear classical and quantum dynam- 
ical systems. This implies to tackle the problem of the 
statistical behaviour of infinite systems. 

A distinctive feature of the radiation matter interac- 
tion in the black body is that each normal mode, or 
field's oscillator, interacts with the matter's degrees of 
freedom only, and that the strength of this interaction 
decreases as the mode's frequency increases. In this pa- 
per we introduce a simplified model which however shares 
the main features of the black body. We consider a sys- 
tem of N oscillators with mass m$ = c/i 2 and frequency 
LUi = yjk/rrii — ai with a — Wkjc. When at its central 
position x — each oscillator collides elastically with a 
particle of mass M » rrii. Therefore the whole system 
is conservative with total energy 
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E = E +J2 E i= (1/2)MV 2 + Im + h. 
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where Eq is the energy of the heavy particle and Ii 
are the actions of the oscillators. The interaction be- 
tween oscillators and the heavy particle is provided by 
the hardcore(h.c) collisions. In our numerical compu- 
tations we have taken the mass of the heavy particle 
M = (V5 + l)/2, c = 0.51, k = 0.1 so that = ai 
with a w 0.443. Notice that our system is of a billiard 
type and therefore the trajectory does not depend on the 
total energy E which is merely a scale factor. We recall 
that the one-dimensional system of two hard core point 
particles with fixed boundary conditions is equivalent to 
the billiard in a right triangle which is assumed to be 
ergodic and weakly mixing jjj. In our model the parti- 
cle exchanges energy with each oscillator and moreover, 
after each collision, we assign at random, the sign of the 
velocity of the heavy particle. Our classical model has 
therefore a high degree of chaoticity. We would like to 
stress that for the purpose of the present paper we only 
need a mechanisms which conserves the total energy and 
allows energy exchange between the different oscillators. 
Therefore, we are not interested in the detailed motion 
of the heavy particle. Instead we need a sufficient degree 
of chaotic behaviour in order to ensure ergodicity of the 
entire system with N + 1 degrees of freedom g] . 

The model system (||) is a mechanical version of the 
one dimensional black body problem discussed in || and 
one expects a similar statistical behaviour. Indeed we 
have numerically computed the time averaged energies 
of the oscillators and of the particle. As expected, we 
have observed that for any N, and independently of the 
initial condition, the system approaches the equipartition 
state in which the time averaged energies of the oscilla- 
tors and of the particle are all equal (an example is given 
by the open squares in fig. 1) . Therefore, for any finite 
value of the total energy E, with increasing the number 
N of oscillators the temperature of the system T = E/N 
will decrease down to zero as N — *■ oo . The mechanisms 
through which such state is approached is the one already 
envisaged by Jeans 0]: energy flows from the matter (our 
particle) to higher and higher modes of the electromag- 
netic field(our oscillators) in such a way that the time 
averaged energy in each mode is zero while the total en- 
ergy remains constant. Therefore the field continuously 
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absorbs energy from the matter, energy will move endless 
to higher and higher modes, and the whole system cools 
down to absolute zero temperature. This would be the 
behaviour according to classical laws. Fortunately how- 
ever, this is not the case since our world is governed by 
quantum mechanics. The latter leads to the celebrated 
Planck distribution which, in one dimension, reads: 

Ek = E{u>i) = hui/(expj3hui - 1) (2) 

The temperature T = 1/(3 is finite and connected to 
the total energy Ef in the field (in the oscillators in our 
case) by the relation Ef — aT 2 which is known as Stefan- 
Boltzmann law (in the one dimensional case). 

Expression^) has been derived by statistical methods. 
A dynamical derivation based on the quantum theory or 
on the numerical solution of the Schrodinger equation 
without any statistical assumption is still lacking. In this 
sense we are asking, for the quantum case, the same ques- 
tion that was posed about 50 years ago by Fermi, Pasta 
and Ulam || regarding the problem of classical equipar- 
tition (which is the limit of @ for % — ► 0). Equiparti- 
tion was a very well established consequence of classical 
statistical mechanics even though it rested on statisti- 
cal assumptions. Fermi however, wanted to derive this 
property by a direct solution of Newton equations. The 
result was unexpected and Fermi considered it as one of 
the most important of his life. Indeed, as we now know, 
it is at the root of the modern field of nonlinear clas- 
sical dynamics and chaos. It is therefore worthwhile to 
ask now a similar question: can we, by direct numerical 
integration of the time dependent Schrodinger equation, 
derive expression^)? 

A convenient guide in this direction can be provided by 
recent progress in the so-called field of quantum chaos. 
One of the main results in the study of this field is the 
quantum suppression of classically chaotic diffusion . 
This phenomenon is known as quantum dynamical local- 
ization since is the dynamical analog of Anderson local- 
ization which takes place in disordered solids ||. Local- 
ization represents a strong deviation from quantum er- 
godicity and only when localization length is larger than 
the sample size then eigenfunctions are extended, quan- 
tum ergodicity takes place and statistical methods can 
be applied. For example, in such situation, statistical 
properties of energy levels can be described by random 
matrix theory. Therefore, if the eigenfunctions are local- 
ized, then quantum evolution leads to a stationary state 
which is typically exponentially localized around the ini- 
tially excited state and therefore strongly depends on the 
initial condition. On the other hand, if eigenfunctions are 
delocalised then the quantum stationary state is close to 
the corresponding classical one. In both cases it is not 
clear according to which quantum dynamical mechanisms 
Planck distribution sets in M . 

In order to provide a clear answer to the above question 
it would be necessary to compute the quantum evolution 



for a system with many interacting particles like the sys- 
tem (|j). This is a too difficult task. 
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FIG. 1. Time averaged energy of the particle Eo and of 
the oscillators Ei , as a function of frequency U) for the discrete 
model with TV = 64 and total energy E — 60. The averages 
are taken over a large number t of total collisions. Initial con- 
dition with almost all energy in the heavy particle: t — 10 7 
(+), t — 4(10 7 ) (*);classical model, t = 4(10 7 ) (open squares). 
Initial condition with almost all energy in the oscillator i=32, 
t — 10 7 (x). The full curve is the Planck law (2) with /3 given 
by eq(3) at E = 60(/3 = 0.252). 

In the following we present instead the results of a 
numerical integration of classical system (|l|) in which 
however we allow for the actions Ii to take only inte- 
gers values. Namely, after each collision, which obeys 
the usual classical conservation laws, we substitute the 
values Ij = Ej/ujj with the nearby integers rij. The 
choice between the upper or the lower nearest integer is 
made at random after each collision. The roundoff en- 
ergy of the oscillator is then given to the heavy particle 
in order to ensure total energy conservation. We will re- 
fer to this model as to the discrete model to distinguish 
it from the usual model(|l|) in which classical evolution is 
computed without any approximation. Our surmise, to 
be verified, is that such a simple discretization procedure 
will qualitatively reproduce the main results of an exact 
integration of Schrodinger equation (with h = 1). Such 
possibility was also suggested in |ic]-|l2]]. More recently 
we have found that application of such procedure to the 
standard map leads to exponential quantum localization 
p3| . Clearly the difference between quantum and classi- 
cal mechanics goes much beyond the discrete nature of 
phase space. It is however our hope to gain a better 
understanding of their relationship. 

In fig.l we show the time averaged energies of the par- 
ticle and of the oscillators, obtained with the above de- 
scribed numerical scheme for N — 64 and total energy 
E = 60. The thin curve is the theoretical Planck law (0) 
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with Ti = 1 and (3 given by (|hwith E = 60. 

oo 

B = l//3 + 5^Wi/(exp/9wi-l) (3) 

i=l 

It is quite remarkable that, independently on the ini- 
tial condition, the discrete model always reaches the same 
stationary distribution. We have also checked that the 
time-averaged values do not change by increasing the in- 
tegration time. As an example we show in fig. 2 the 
time-averaged energies for few oscillators as a function 
of time measured in number of collisions, up to t = 10 7 . 
These values do not change by further increasing the in- 
tegration time. 

We will not enter here in several questions of detail and 
simply notice that the overall agreement between the nu- 
merical results and the theoretical curve (Q) (with the 
correct value of fi — l)is sufficiently satisfactory. Solu- 
tion of the classical model starting with the same initial 
conditions leads to equipartition as is shown by the open 
squares in fig.l. 
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FIG. 2. Time-averaged energy < Ei > of oscillators as a 
function of total number t of collisions. The curves refer to 
oscillators (starting from above): i = 1, i = 8, i = 16, i = 24, 
i = 32, i = 40. 

In order to analyse the behaviour of the system as a 
function of the energy, at fixed N, we plot in fig. 3 the 
total energy E of the system over the time averaged 
energy< Eo > of the particle. Clearly in our discret 
model the classical limit is reached for E — > oo, since the 
effects of discretization become less and less important as 
the total energy is increased. Therefore, for sufficiently 
large energy, one expects equipartition which is given by 
the relation < E >= (N + 1) < E >. According to stan- 
dard statistical mechanics, the temperature is defined as 
the average kinetic energy T =< Eq) >. 
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FIG. 3. The total energy E of the system over the aver- 
age energy < Eo >= T of the particle for different TV values. 
N — 64(open rombus), TV = 32(+), N = 16(open squares), 
N = 8(x). The solid line has slope two while the four parallel 
lines have slope one. 



As we decrease energy, at fixed particle number N, 
we expect deviation from equipartition. It is really re- 
markable that integration of our discrete model leads to 
the Stefan-Bolzmann law which, in our one-dimensional 
case, implies that the total energy increases quadrati- 
cally with the temperature. The solid line in fig. 3 is 
given by E = o~< E > 2 = aT 2 with a w 6 (indepen- 
dent on N). The four parallel lines are given instead by 
E = (N + 1) < E > {N = 8,16,32,64)which is the 
exact equipartition law of classical mechanics. Accord- 
ing to our data, the transition from Stefan-Boltzmann to 
equipartition is quite sharp and takes place at E » N 2 /a. 
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FIG. 4. Plot of E/ < Eo > versus TV for different val- 
ues of total energy E; E = 25(open triangles), E — 100(x), 
E = 225(open squares), E = 400(+), E = 1600(open rom- 
bus). The dotted line has slope one and corresponds to 
equipartition. 
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In order to analyse instead the behaviour of the system 
as a function of the number N of oscillators, at fixed total 
energy E, we plot in fig. 4 the quantity Ej < E n > over 
N, for different values of the total energy E. The dotted 
line corresponds to equipartition and it is approached by 
decreasing TV, at fixed energy E. Instead, by increas- 
ing N one approaches a constant value oc \f~E~. This is 
again in agreement with the Stefan-Boltzmann law. In 
summary, the results presented in figs 1-4 agree with the 
predictions of quantum mechanics. For a given N, at low 
temperature, the Stefan-Boltzmann law and the Planck 
distribution are obtained; as temperature is increased, 
one approaches the semiclassical limit in which equipar- 
tition takes place. 

It may be interesting to examine this transition in re- 
lation to the general theory of quantum dynamical locali- 
sation. According to [H, the scaling Ansatz is equivalent 
to postulating the existence of a function f(x) such that 

l/N = f(x),x = t/N (4) 

where I is the actual localization length in the sam- 
ple of finite size AT, and £ is the localization length for 
the infinite sample. As a measure of localization, follow- 
ing standard procedure, we take the inverse participation 
ratio: 

N 

l=<E f > 2 /J2<E j > 2 (5) 
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FIG. 5. The actual localization length l/N for the dis- 
crete model over the localization length £/JV for the infi- 
nite "sample" for different "sample sizes" N and for dif- 
ferent energies E. N = 64(+), N = 56(x), N = 48(*), 
N = 40(open squares), N = 32(full squares), TV = 24(open 
circles), N = 16(full circles), N = 8(open triangles). 

In fig. 5 we plot the actual localization length l/N ver- 
sus £/N where £ is computed from the distribution (||). 



The localisation length I is obtained from expression (|5J) 
where < Ei > are the time-averaged energies of the os- 
cillators computed for the discrete model and < Ej > is 
the total average energy of the oscillators. A quite nice 
scaling behaviour is observed. 

Certainly the results presented here require a better 
understanding. In particular it is not clear to what ex- 
tent they are related to the peculiarity of the radiation 
matter interaction in which the field modes interact only 
through the matter. Needless to say the conclusions we 
have drawn here refer to the discrete model ([!]). It is our 
personal opinion that the solution of the time dependent 
Schrodinger equation will lead to the same qualitative re- 
sults of figs. 1-5. This however needs a carefull study of 
true quantum mechanics of system (|l|) or of similar mod- 
els. 
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